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FIRST ORDER STEP RESPONSE IDENTIFICATION FROM 
NOISY DATA 


Dan STEFANOIU!, Nicolai CHRISTOV’, Vasilica VOINEA! 


Abstract. In the real world, systems and signals are affected by stochastic perturbations 
briefly referred to as ,,noises”’. Such noises can be generated by unknown sources located 
either inside or outside of a system and usually corrupt in unknown way the acquired 
data the system can provide. Depending on the noises power in the acquired data, some 
characteristics of the system under study can or cannot be determined. This article 
introduces a method to identify optimal smooth step response, of first order, from noisy 
data, by means of Newton-Raphson method employed to minimize a quadratic criterion. 
Simulations with real world data prove the method effectiveness. 
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1. Introduction and problem statement 


Automatic Control is a scientific and engineering field that strongly relies on 
Systems Theory [7], [3], [6], System Identification [8], [4] and Optimization 
Theory [1], [9], as main pillars. Since its inception, despite the very complex 
theoretical approach that the above-mentioned theories have reached, a 
humongous number of automatic control applications were developed. Usual 
practitioners are mostly interested in applying rather simple theoretical results to 
real world systems. Nevertheless, lately, in some applications, there is an 
increasing interest of filling the gap between advanced (often complex) theoretical 
results and numerical procedures that can be implemented on their basis. For 
example, in Industrial Automation field [5], very seldom controllers outside the 
PID class are accepted, whilst, in Aerospace Industry [2], optimal state space 
controllers are already implemented, although their complexity is sensibly higher. 

This article tries to answer a question of interest for many Automatic Control 
practitioners, namely: how to extract some characteristics of a dynamic system 
from real world, by using acquired data the system can provide? At a first sight, 
the problem related to this question should not be so difficult to solve. 
Nevertheless, complications arise because the signals of real world are corrupted 
by stochastic perturbations, also referred to as noises. Unfortunately, in most 
cases, one cannot know what the noise sources are, whether they are located 
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inside or outside the system under study and how the noises are mixed with the 
signals to acquire. Or, especially in Systems Theory, many theoretical results are 
concerned with dynamical entities that evolve in a clean environment, either noise 
free or with very small influence of noises. Albert Einstein expressed long time 
ago a sad thought: “It always distresses me to see how a beautiful theory is 
destroyed by an ugly reality...”. Hopefully, although not perfect, there are 
scientific tools allowing the user to draw some line between the useful 
information the acquired data encode and the parasite information induced by 
noises in the data. 

In this article, only a small part of the general problem stated above is 
approached. Expressis verbis, one seeks to build a smooth first order step response 
from a dataset corrupted by noises. Such a response is displayed in Figure 1 and 
exhibits the variation of air flow temperature acquired at the output of a 
climatizer. As one can easily notice, this response suggests that, maybe, the 
climatizer can be modelled as a first order system that responded to some step 
input signal. However, the theoretical (smooth) first order step response is not 
easy to draw empirically. 
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Figure 1. Noisy temperature data acquired from a climatizer output 


The smooth step response could help the user to extract some practical features 
of the dynamical system that yielded acquiring the dataset. More specifically, 
from such a response, one can estimate the intrinsic delay T, the time constant 7’, 
the DC gain K and the initial offset Ay, as shown in Figure 2. Subsequently, the 


4 parameters can help the user to tune an appropriate or optimal controller for the 
climatizer (e.g. of PID type, by using Broida method [10]). 

To determine the best step response for a given dataset, an optimization 
criterion is needed. In this case, the usual quadratic criterion can be employed. Its 
minimization can be performed by means of various exact optimization 
techniques [1]. 
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Figure 2. Theoretical first order step response and its main parameters 


Such a technique is the Newton-Raphson’s (N-R) one, due to the possibility to 
compute not only the gradient (or the first derivative) of quadratic criterion, but 
also its Hessian matrix (or the second derivative). 

The article is structured as follows. The next section is devoted to 
mathematical formulation of optimization problems to solve. In section 3, one 
shows how the optimal solution can theoretically be obtained, by means of N-R 
numerical procedure. Some simulation results on real world data are presented 
and analyzed in Section 4. The article completes with concluding remarks and a 
references list. 


2. Optimization problems to solve 
One starts with the mathematical expression of smooth (noise free) first order 
step response, in continuous time: 


Ay »t €[0,T) 


(t)= =A : (1) 
z K(1-e "| .Ay $t2T 
where the shape is uniquely determined by 4 parameters: t>O (the intrinsic 
delay), 7 >O (the time constant), K €IR’ (the DC gain) and Aye R (the 
initial offset). This shape should fit to the best to the available dataset 
dD = ees with N EN’. Obviously, the data were acquired with some 


sampling period, T >Q (unit, if unknown). Consequently, before anything else, 
the variable part of step response (1) has to be discretized with the same period: 


nT, -t 


wnl=y(nT)=K(1-e 7) | +ay, Vn>n = A (2) 


s 
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Then, the following quadratic criterion can be employed to determine the best 
fitting step response to the dataset dD, : 


YU, (,T, K,Ay) => '(y,-Ay) + Sy, = K(I-eF }-ay| » 16) 
It follows that the optimization problem can be formulated as follows: 
(4.7,K.Ay)= argmin %,(t,T,K,Ay). (4) 


In general, the solution of problem (4) can be found by solving the null 
gradient equation: 


VV,(t,T,K,Ay) =0. (5) 
However, the minimum of VU. is obtained when both sums of quadratic errors are 


minimum. This means the first sum is minimum when: 


n,-l 


d [E(v,-ayy f=0 ee -25(y,-ay)=0. ‘C\6 


o(Ay) n=0 n=0 
Thus, the optimum initial offset is: 
| nz-l 
AY == (7) 
nN, n=0 
once the optimum intrinsic delay T has been found. This corresponds well to the 
image in Figure 2. 


The other optimum parameters are solutions of the compatible system below, 
obtained from (5): 


’ ely, =x(1Se'e" ay) =0 
Si(1-e'e*] y-K( 1-6” )-ay =o. (8) 


Since the system (8) is nonlinear, the optimal solution is not easy to compute. 
Therefore, it is suitable to simplify the approach and to search for a sub-optimal 
solution, with acceptable accuracy. 


This allows splitting the optimization problem in two sub-problems, by 
grouping the unknown parameters in two couples: {K ,Ay} on one side and 


{t,T} on the other side. Two new quadratic criteria can now be defined. The 
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equations of system (8) suggest that, to simplify the solving procedure of 
optimization problems, the following notations should be employed: 


E=e €(0,1) and Q=e D1. (9) 


Then, the new optimization criteria are: 


4,(K) => (y,—ay)-K(1-08)] : (10) 
4;(0.8)= 3] (9, -ay-K) + ROE]: ay 


In definition (10), one considers that the couple {2,7} is already estimated, 


which allows computing 6 and E according to notations (9). Also, the discrete 
time delay is n,=[%/T. ]. Similarly, in definition (11), the DC gain K is 
known. In both criteria, the offset Ay is computed by means of equation (7). In 


criterion (11), the delay T will be estimated recursively, so that Ay can always 


be.computed from the delay value available at previous iteration. 
The criterion (10) can equivalently be expressed by performing the following 
index changing: n <~n—n,. Thus, 
N-n;-1 Kx 2 
U.(K)= > Moe ay) 2K (1 6" )| ; (12) 
n=0 


Moreover, equation (12) reveals that, instead of using the truncated data 


{Youn \ oe? One can work with the modified data set: 


Ven = Varn ~AY, Vne0,N, -1, (13) 
where NV, = N —n.. Then, the criterion (10) becomes: 
4.(K) =>] y..-K (1-68) | 


which simplifies its definition. The new definition, (14), involves that two 
preliminary operations have to be applied on the genuine dataset: (a) remove the 
subset located into the delay zone delimited by n,; (b) shift the remaining data by 


2 
’ 


(14) 


the offset Ay, as illustrated in Figure 3. 
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Figure 3. Preliminary operations to apply on dataset for the quadratic criterion (10) 


In case of criterion (11), the approach is similar. By performing the index 
changing n <— N—/N_, a new expression 1s obtained: 


V, (8, S) = > ( eo Ay - k) + Keg" | : (15) 
According to (15), the modified data set is: 
Wn = Nie = ALF K , Vne 0,N.-1, (16) 


where NV —— N- Nn. is the variable number of data to work with. Practically, the 
operations to apply on the genuine dataset are almost the same as pointed in 


Figure 3. One difference is that, here, the data are shifted by the sum Ay + K 


(instead of Ay ), which enforces the steady-state part to vary around the time axis 


(provided the offset and the DC gain were correctly estimated). Another 
difference is given by the delay Tt, which, here, is unknown. 
From equation (15), one obtains: 


4; (0,2)=>(9..,+ ROE") = 
7 oo Mel "7 (17) 
=(Feno + ROE") + Y(F.c4 + ROE). 


n=l 


The final expression in (17) allows making the difference between two cases: 
n.=0 and n_ > 0. If no intrinsic delay exists (t = 0), then the first term of final 


expression in (17) is independent on E (as .. =1); also, the criterion VU: does 
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0 


not depend on 0 variable (as 8 =e’ =1). Otherwise (t>0), the full sum in 


(17) can be employed (as 0< &" <1 and 0=e’ >1). 
The optimization problems to solve are then: 


K =argmin@'(K); (18) 
(6,2) =<argmin % (0,8). (19) 
@>1,E€(0,1) 


3. Solving the optimization problems 


Like in case of overall criterion (3), to solve the optimization problems (18) 
and (19), the gradients of criteria (14) and (17) need to be computed in the first 
place. Then, the solutions result by solving the equations obtained when enforcing 
the gradients to null values. 


After some simple algebraic manipulations, the equations to solve are: 
N;-1 


> (1- 62" )| ,, K-68 )| =0. (20) 


Sees (5... + KOE") =O" 50> 6>1). (21) 
dine (y, + KB) =0 LAND (t=0> 0=1) 
- (22) 


N,-1 


(mtn end. + Ko") =0 ,n.>0 (t>0>051) 


n=0 
Note that, in no delay case, since the criterion (17) only depends on € variable, 


the system (21)-(22) is replaced by the first equation of (22). 
Focus first on equation (20), which directly gives: 


K==. : (23) 


To yield efficient implementation of solution (23), the sums of geometric series in 
denominator can be computed in advance. Thus: 


N;-1 ra 


5 (Ge) =n, 26888 ee 


n=0 =0 


sek © N; xt, six 2Ne 
268" S "5 Gre S_—* 
ool Gc 


I] 
= 
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2, 6 If ac, EY +1 
=N,+ 6é" a i [ Gc" = } 


(24) 
E+1 
Approach now the system (21)-(22). 


a. If t=O (0=1), only the first equation of (22) is available 


(25) 


The sum in the right side of equation (25) can be computed by using 
geometric series. Hence: 


N-1 


ne tile 89S) S765 Pa 


7 “a(E) ef 
(26) 
(N -1)&" — NG"? +4 


By inserting (26) in (25), the equation to solve becomes: 
N-1 hy 
(6-1) Dinky, =-KE| (N -1)E" = NE"? +1]. (27) 
n=l 
b.If t>O, the product KO can be expressed from (21) and second equation 
of (22): 


) 
co (Eytan Le Ses... (Ze 
eo of Sng Se... -(Se)8 nh 


(28) 
n=0 n=l 


Equation (28) can be transformed by means of geometric series (see 
property (26) as well): 


é[ (NW, -1)e" - NE +1] DES Von = 


n=0 
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N,- 


=(€°-1)(6" -1) SB", 29) 


n=l 
None of equations (27) or (29) can be solved in closed form, as the unknown 
variable € cannot be extracted from the sums with terms including data samples. 


Therefore, the N-R procedure can be employed to approximate both solutions. 
Recall that the N-R procedure operates with gradients and Hessian matrices of 

a cost function. In case of equations (27) or (29), the corresponding cost function 

only depends on scalar €. The first two derivatives of cost function are sufficient 


to run the procedure, without actually knowing the cost function itself. 
Fortunately, each equation leads to the definition of cost function first derivative. 
More specifically: 


a.if t=O, then: 
FG) =(8 1) Sombery, FRE[(N —1)B" — NE" +1]; GO) 


b. otherwise: 


N,-1 


fO=(-Y)(e"-)py nega 
a (31) 


~é[(N,-DE* -NEN HIS es... 


n=0 


for any T>O. 
From (30) and (31), the second derivatives can be computed straightforwardly: 


a.if t=O, then: 


Fete) = 28 (e? - ne" y+ (e*=1)° Fann Ne" yt 


} (32) 
" K[(W ~1)(2N +e" —N(QN De"? + 1]; 


b. otherwise: 


FXG) =E[(N,+3)E" NE" —3]S ENS. + 


+(e 1G" $n DEB... 
—& 2" L(N. —h(2n, a aN (EN, -H X85... o 
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Using the N-R procedure is quite easy in case of null intrinsic delay. However, 
if tT>QO, the main problem is the dependency of both derivatives on that delay. 
Therefore, to allow implementation of N-R procedure, the following strategy 
could be adopted. 


A. If the initial time delay is null, then only equations (30), (32) and (23)-(24) 
can be employed. In this case, n_ always is null and only parameters 
{K,T} can be estimated. 


B. If the initial time delay is non null: 


1. Use the (initial) value of t > O to find the (first) value E, by means of 
N-R procedure (equations (31) and (33)). 


2. With &, the parameter @ can be estimated by using equation (21): 


6= 2 ae 2 a (34) 


3.From é and 6 , the delay can be estimated with equations (9): 


jel 3 t=T log6. (35) 
logé 


4. Estimate the initial offset Ay by means of equation (7). 


5.Use equation (23) to estimate the DC gain, K : 
6. Repeat steps 1-5 with until the accuracy condition is met. 


The model performance can be assessed by means of criterion (3) or, even 
better, by means of signal-to-noise ratio (SNR), computed with the help of 
criterion (3): 

1 ales ay) 
er ee) 
SNR (#7, K,Ay)= “ve (36) 
4. #,7,KAy) 
N 


dB 


1 es 


In definition (36), y is the noisy data average. Also, the SNR is expressed in 


decibels (dB). (Recall that al, = 201log(a) .) Thus, the SNR is the ratio between 


the standard deviation of data and standard deviation of modeling error. The 
higher the SNR, the better the model. 
The algorithm to solve the optimization problems is listed next. 
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Algorithm to find optimal step response, based on Newton-Raphson procedure 


> Inputs: 
eA noisy dataset: D, = { y, _—— ,with NEN’. 
e Sampling period: 7’ > 0 (unit, if unknown). 
e Initial values of parameters: K,, T,, T,. 
e Accuracy threshold: € (equal to 10, by default). 


Initialization 
«Compute the discrete-time delay n, = Ee / T | : 


"Compute &, =exp(-T./T,). 
= Set the initial variable step in N-R procedure: @, =1. 
Alive, =O (=n, =0): 
— set t=O and Ay = 0, as they cannot be estimated, in this case; 
— compute f/(&,) with (30) and f(&,) with (32), for K=K;; 
ACE) 
ACE.) 


— estimate the first correction: A .= 7 
= Otherwise (t, >O =n, > 0): 


— compute NV, = N —n, and Ay = LL 


—Ay,-K,, Vne0,N, -1; 
— compute f'(&,) with (1) and f"(&,) with (3) for N.=WN, and 


A 


YeRn ae own? Vne0,N, —1. 


— extract the data beyond n,: y,, = y 


n+ng 


— estimate the first correction: A, =— 


= Set the initial number of iterations: k =O 

If t, =0, do: 

1. Approach the optimum: €,,, =&, + a,A,. 

2. Update the DC gain K,.. with (23)-(24), for €= 
case, 9=1 and N.=N.) 


. (Note that, in this 


k+l 
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3. Compute f/(&,,,) with (30) and f,(&,,,) with (32). (Note that, in this 
case, K= K,.,.) 
AiG) 


5. Update the variable step: O,,, =, + —* 


k 


4. Update the correction: A, =— 


k 
6. Increment the number of iterations: k <—k +1. 
until \A,| <Se 


If t, > 0, do: 


1. Approach the optimum: 6, =€, +a,A,. 


with (34); for E=€ 
n.=n,and J.,, Fy) Vne0,™% -1) 

3. Update JT. =—-T /logé,,, and t,,, =7,,, log® 

. Update: n,,, = co /T | and NV. =N-n.,,. 


1 Mei =l 
. Update the initial offset: Ay,,, =—— >. y,. 
nN 


n=0 
k+l 


na 


. (Note that, in this case, K = K_, 


k+l k 


2. Compute 0,, 


1 


k+1* 


. Extract the data beyond n,_: 
Vilin = Jaen, To AV er 7 Vir €0,N,,, - 1. 
. Update the DC gain K,.. with (23)-(24), for €=€ 
N.=N.,,,- 
. Compute f’(&,,,) with 31) and f"(&,,,) with (33). (Note that, in this 
case, N.=N. =Vore Vne 0,N,-1) 
FACS 
Ge 


A 
10. Update the variable step: &,, =O, + a ; 


6=60., and 


k+1? 


1 and Ves 


. Update the correction: A, =— 


k 
k 


11. Increment the number of iterations: k <k +1. 
until A, | KE. 
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“+ Final computations 
"Set: K=K,, T =T,,. 
If t, >0,set: T=1,, Ay =Ay,. 
“Compute the optimal step response ) with (1) and (2), by using the 
sampling period T beyond [4 J Eats 


= Estimate the step response performance by computing SNR Ge K : Ay) 


with (36). 
Outputs: 


na aA 


e Optimal values of parameters: K , 7, t, Ay. 


eOptimal step response, sampled with T : Dd. =) S ___, with same 


ne0,N-1” 


length as the initial dataset ( JN ). 


e Step response performance: SNR(#7, K Ay) (in dB). 


e Number of iterations to meet the accuracy condition: k . 


4. Simulation results and discussion 


As already mentioned in Introduction, the dataset to test was acquired at the 
output of a climatizer mounted in a large room, with many sources of noises (such 
as doors and windows that are opened and closed very frequently). The 
temperature of airflow was measured and recorded with the sampling period 
T =1 min. The number of samples is N = 10041 (quite big). Figure 1 displays 


the variation of temperature. One can notice the noises are quite powerful, so that 
it is difficult to manually draw a fitted step response. 

The algorithm above has been implemented into MATLAB™ programming 
environment. The procedure was initiated to run with the initialization below: 


K,=1°C, T, =430 min, t, =1000 min, ¢=10", (37) 


which has been empirically obtained, like in Figure 3, after visually inspecting the 
variation in Figure 1. The optimization results are illustrated in Figure 2, where 
the optimal parameters are: 


K =1.043°C, T = 404.5 min, t= 1367.7 min, Ay =21.89°C. (38) 
Moreover, the steady-state temperature can also be estimated: 


K + Ay = 22.93°C (39) 
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Apparently, the climatizer tries to keep the room temperature as close as possible 
to the value of 23°C. 

The performance of N-R procedure is revealed not only by the small number of 
iterations until the optimal values were found, k =12, but also by the search 
duration, AJ =0.32s (as obtained on a personal computer from i7 family, of 


10-th generation, with 8 parallel processors). 
The performance of optimal step response is given by: 


SNR( 27K, Ay) ~12,987dB. (40) 


The main limitation of the algorithm above is its sensitivity to initialization. 
The simulations revealed that the procedure is the most sensitive to small 
variations of initial DC gain values. For example, if the initial DC gain decreases 


to K,=0.9°C (with 0.1°C only), while keeping the other initial values, the 


sub-optimal result of Figure 4 is obtained. Comparing to the optimal result in 
Figure 2, this step response has smaller performance, which can be observed 
easily. 
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Figure 4. Sub-optimal step response, obtained after decreasing the initial DC gain by 0.1°C. 


Moreover, the SNR decrease confirms the visual observation: 
SNR(#,7,K,Ay)=7.NdB. (41) 


The performance deteriorated with more than 45%, although the initial value of 
DC gain decreased by 10% only. So, it is important to correctly set the initial 
values of parameters. To this respect, as a future development, one can use a 
metaheuristic approach [9] to solve the following optimization problem: what 
initialization leads to the highest SNR? 

When looking back at the variation of temperature, even more closely, one 
realizes that, beside the first jump of about 1°C, in the steady-state zone, the 
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internal regulator of climatizer seems to be a rudimentary one, working in steps. 
Thus, the temperature variation includes several small successive step responses, 
as reaction to that regulator. Then, the optimization algorithm can be employed to 
build the noise free response of climatizer, under the control of regulator. 
Evidently, this attempt involves building a set of local initializations, one for each 
zone where a small step response seems to exist. Although this activity could be 
tedious, if the initializations are realistic, the reward is guaranteed. 

After running the algorithm for each possible local step response, the result is 
displayed in Figure 5. 
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Figure 5. Optimal step responses, obtained after running the N-R procedure on each local variation 
zone of temperature. 


The smooth variation shows that the climatizer has different parameters of step 
response, depending on the height of reference step (or the setpoint) to follow. 
This behavior was expected, as the static characteristic of such a climatizer is 
nonlinear and exhibits hysteresis phenomenon. Hence, changing from a nominal 
point to another, on the static characteristic, automatically involves changing the 
parameters of step response (especially of time constant). For the multiple step 
responses in Figure 5, denoted by. , the performance is: 


SNR (5) = 20.54dB. (42) 
There is no wonder the SNR has sensibly increased, comparing to single step 


response case, as the multiple step responses curve better fits the dataset. 


Concluding remarks 


This article approached the problem of optimal first order step responses 
estimation by using noisy datasets. The proposed method to solve the problem 
relies on quadratic criteria and Newton-Raphson procedure. The simulation results 
have proven both the method effectiveness and its main limitation, namely the 
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high sensitivity to initialization of some parameters (especially of DC gain). This 
drawback suggests using a metaheuristic to find the initialization that maximizes 
the SNR. Another possible future development is to design a method allowing the 
user to find the optimal second order step response from noisy data. It is well 
known that smooth first and second order step responses are extremely useful in 
Systems Theory and Industrial Automation. 
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